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Abstract 

We present a sophisticated gamma-ray likelihood reconstruction technique for 
Imaging Atmospheric Cerenkov Telescopes. The technique is based on the compar- 
ison of the raw Cherenkov camera pixel images of a photon induced atmospheric 
particle shower with the predictions from a semi-analytical model. The approach 
was initiated by the CAT experiment in the 1990's, and has been further developed 
by a new fit algorithm based on a log-likelihood minimisation using all pixels in the 
camera, a precise treatment of night sky background noise, the use of stereoscopy 
and the introduction of first interaction depth as parameter of the model. 

The reconstruction technique provides a more precise direction and energy re- 
construction of the photon induced shower compared to other techniques in use, 
together with a better gamma efficiency, especially at low energies, as well as an 
improved background rejection. For data taken with the H.E.S.S. experiment, the 
reconstruction technique yielded a factor of ~2 better sensitivity compared to the 
H.E.S.S. standard reconstruction techniques based on second moments of the cam- 
era images (Hillas Parameter technique). 
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Introduction 



The last decade saw the emergence of very high energy (VHE; E > 100 GeV) 
gamma-ray astronomy as a new observational discipline, with the number of 
VHE gamma-ray sources now approaching 100. This success was driven by the 
third generation of ground based Imaging Atmospheric Chcrcnkov Telescope 
Systems (lACT), such as H.E.S.S., with an order of magnitude better sen- 
sitivity compared to the previous generation instruments (see e.g. [1,2] for a 
recent review). The improvement was made possible by the use of stereoscopic 
systems of large telescopes equipped with finely pixelated fast cameras. 

To reconstruct the direction and energy of the primary gamma-ray and to dis- 
criminate them from charged cosmic rays most of the current experiments use 
reconstruction techniques based on the second moments of the pixel ampli- 
tudes in the camera (Hillas parameters) [3,4]. These techniques are very robust 
and efficient. However, a more sophisticated albeit more computing time in- 
tensive reconstruction technique was pioneered by the CAT experiment [5], 
taking advantage of its very finely pixelized camera. The technique is based 
on a comparison of the recorded Cherenkov light distributions of a photon 
induced electromagnetic shower in the camera, i.e. the shower images, with 
calculated shower images from a model of the Cherenkov light distribution in 
electromagnetic showers. The reconstruction technique presented in this pa- 
per [Model Analysis) is a further development and improvement of this early 
work. 

The calculated shower images are derived from the Cherenkov light distribu- 
tion of charged particles in electromagnetic showers taking into account light 
collection efficiencies, atmospheric absorption etc. The Cherenkov light dis- 
tribution of a shower is determined by the longitudinal, lateral, and angular 
distributions of charged particles in the shower. These distributions are de- 
rived from Monte Carlo simulations and parametrised to yield an analytical 
description of the shower, i.e. the shower model, including the depth of the 
first interaction as a new parameter in the parametrisation. Additionally, the 
contribution of the night sky background noise in the camera in every pixel is 
modelled on the basis of a detailed statistical analysis. Thus the fit procedure 
does not require a dedicated image cleaning procedure to extract the pixels 
illuminated by the shower. The parameters of the calculated shower that best 
fits the measured shower image are determined in a minimisation procedure 
which yields a selection criteria to discriminate the gamma-ray induced shower 
from the hadronic background. 

The different parts of the model, i.e. the semi-analytical description of the 
shower development are described in section 1, and are used in section 2 to 
generate the shower image templates. In section 3 the fit algorithm is de- 
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scribed. Finally, section 4 is devoted to the description of the performance 
of source analysis using this model reconstruction, with detailed comparison 
with alternate reconstruction techniques. 



1 Charged particle distributions in an electromagnetic shower 



The model of the pixel illumination by Chcrcnkov light in the camera is based 
on analytical descriptions of the energy-dependent longitudinal, lateral, and 
angular distributions of charged particles in electromagnetic showers. The elec- 
tromagnetic showers were generated with the KASKADE [6] shower simula- 
tion program. The KASKADE program has been improved to include among 
others a more precise treatment of Bhabha and Moller scattering [7,8]. A 
detailed comparison of several generators in the H.E.S.S. -collaboration, in- 
cluding CORSIKA [9], showed no noticeable difference between the different 
generators. 

In order to cover the dynamical range of the H.E.S.S. instruments, showers 
were generated at energies of 10 GeV, 50 GeV, 100 GeV, 500 GeV, 1 TeV, 
5 TeV, 10 TeV and 20 TeV. Typical number of showers required for the gener- 
ation of a smooth model range from ~ 200000 at low energy to a few hundred 
at the highest energies. The varying number of showers ensure that different 
primary energies have similar weights in the model construction. 

The model for atmospheric electromagnetic showers presented here is an ex- 
tension of the model proposed by Hillas in [10]. The nomenclature used here 
follows closely the nomenclature of the model by Hillas. 

As the dominant source of shower-to-shower fluctuations, the depth of the first 
interaction is included as a new parameter in the parametrisation. This helps 
reducing the discrepancy between actual shower images and model predictions, 
and thus improves the reconstruction and discrimination efficiency. 



1.1 Longitudinal distribution of charged particles 



The average number of charged particles J\fe{y, t) in simulated electromagnetic 
showers as a function of the distance to the first interaction point t for different 
primary photon energies are shown in fig. 1. The longitudinal distributions are 
modelled by a modified Greisen formula[ll]: 



K{y,t) = — X exp 
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where y — \n{Eprim/ Ecru) is the scaled primary photon energy with Eprim 
being the primary photon energy and E^-ru being the critical energy {Ecru ~ 
83 MeV in air). The first part of the expression corresponds to the standard 
Greisen formula, with parameters a and b left free. The second part describes 
the decay of the two charged particles created at the first interaction point 
t = and ensures J\fe{y,t = 0) — 2. The shower depth t is given in units 
of radiation lengths. The shower age s measured from the point of the first 
interaction is given by 

, ^ ^ (2) 

1 + cx {b-l)/t ^ ' 



This expression is constructed so that the shower age being s = at the first 
interaction and s = 1 at the shower maximum. The parameter c is the depth 
of shower maximum measured from the first interaction {s — 1 ior t — c) 
and 6 is a scahng factor related to speed of the shower development. The 
parameters a and c are found to be linearly dependent on the scaled energy 
y. A fit of the analytical description on the simulated distributions yielded for 
the parameters a, b and c: 

a = 1.05 + 0.033 x y, 6 = 2.66, c = 0.97 xy- 1.32 (3) 



As can be seen in fig. 1 the analytical description is in good agreement (at a 
level of ~ 5% in the central depth range, degrading to ~ 20% in the shower 
tail and ~ 40% in the shower head) with the simulation for primary energies 
ranging from 10 GeV to more than 20 TeV,with slightly worse performances 
at the very high and very low energies. 

The Cherenkov light distribution of an electromagnetic shower depends on the 
energy dependent longitudinal distributions of charged particles in the shower 
(We/dE with 



K = J ^dE., K{E >Eo) = l ^dE. (4) 

Eo 



Examples of the longitudinal distributions of charged particles dMe/dt{E > 
Eq) in the shower integrated over charge particle energy E above some thresh- 
old ^0 (10 MeV, 20 MeV, 40 MeV, 80 MeV, 150 MeV, 300 MeV, 600 MeV, 
1 GeV and 2 GeV) and for different primary photon energies is shown in 
Fig. 2 . The distributions are modelled using the same analytical expression 
as in Eq. 1 but with a different set of parameters: 



a = (1.058 + 0.014 x y) + 1.6 x 10"^ x |ln£;Mev - 6 



1.5 
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Fig. 1. Top: Longitudinal shower development measured from the first interaction 
point: Number of charged particles in the shower as a function of the distance to the 
first interaction point for several primary energies ranging from 10 GeV to 20 TeV. 
The black histograms show the results of a simulation compared to the analytical 
expression of eq. 1 (solid red curve). Bottom: Ratio of simulation to prediction from 
eq. 1 for various primary energies, showing that eq. 1 is a good description at the 
level of ~ 5% in the central depth range, degrading to ~ 10% at the lowest energies. 

6 = 2.55 + 0.067 X In EmcV (5) 
c = 0.97 X 1.43 + 0.137 x In^MeV - (0.0712 + 0.0005 x y) In^^MeV 



where EueV is the charged particle energy, expressed in MeV. 

The comparison of the analytical function with the distributions from simula- 
tions with different primary energies is shown in Fig. 2. The agreement is very 
good (at the level of ~ 5%) up to particle energies of a few GeV above which 
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Fig. 2. Top and Middle: Energy dependent longitudinal shower development for 
primary energies of (top to bottom and left to right) 10 GeV, 100 GeV, 1 TeV and 
10 TeV, compared to analytical approximation from eq. 5 (solid red lines). Each line 
gives the number of charged particles above some energy (10 MeV, 20 MeV, 40 MeV, 
80 MeV, 150 MeV, 300 MeV, 600 MeV, 1 GeV and 2 GeV from top to bottom) as 
function of atmospheric depth. Bottom: Ratio of simulation to analytical prediction 
from eq. 5 for primary energies of 100 GeV and 10 TeV. 
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the particles contribute very little to the overall Cherenkov light distribution. 



1.2 Angular distribution of particles in a shower 

The angular distribution of charged particles in the shower together with the 
velocity dependent Cherenkov-angle determine the angular distribution of the 
Cherenkov-photons. As described in Moliere theory, the relevant angle of the 
charged particles is the reduced angle w given as 



where 6 is the angle between the direction of the primary particle and the 
direction of the charged particle of energy E in the shower. The angular dis- 
tribution of the charged particles in the shower is decomposed into the mean 
reduced angle (w) as a function of the particle energy of the charged particle 
and the distribution of reduced angles around the mean reduced angle. 

The mean reduced angle is found to depend only very weakly on the primary 
particle energy and is parametrised by 



The mean reduced angle for different primary particle energies as a function 
of the charged particle energy in units of the primary particle energy together 
with the parametrisation is shown in Fig. 3. There is a good overall agreement 
within a few percent between the mean reduced angle from the simulation and 
the parametrisation, in the range of small kinetic energies compared to the 
primary particle energy. Fig. 3 right shows the ratio between the simulation 
and the parametrisation. A deviation is seen when the charged particle kinetic 
energy exceeds about 5% of the primary particle energy (corresponding to a 
red vertical line in fig. 3, right), affecting only a very small fraction of the 
particles in the shower. 

The detail modelling of the Cherenkov light distribution in a shower requires 
the dependency of (w) on the shower age s to be taken into account. The 
parametrisation of eq. 7 is kept, but with parameters varying with shower 






0.435 



(7) 
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Fig. 3. Left: average value of the reduced angle cosine {w) as a function of charged 
particles energy. The solid line is the model parametrisation from equation 7. Right: 
Ratio between simulation distribution and model prediction as function of charged 
particles kinetic energy (scaled to primary particle energy). 

age: 



(w) is) 




P2{s) 



The parameters po, pi and p2 are given by (fig. 4): 



Po{s) = 0.506 X exp(0.351 x In s - 0.147 x In^ s) 

pi(s) = 144 X exp(0.33 x Ins - 0.11 x In^ s) (9) 
P2(s) = 0.887 X exp(0.0507 x Ins - 0.014 x In^ s) 




10 10^ 10^ 10^ 10 10^ 10^ 10^ 10 10° 10^ 10^ 

E [MeV] E [MeV] E [MeV] 



Fig. 4. Mean value of scaled angle (w) (s) as function of particle energy and shower 
age, for different primary particle energies (100 GeV to 10 TeV from left to right). 
The solid lines correspond to the parametrisation of eq. 8. The average scaled par- 
ticle angle is found to be mostly independent of primary particle energy. 



8 



The distribution of angles around the mean is found to be quite independent 
of shower age or primary particle energy. The rescaled angle u is used here: 



The distribution of dlogiQ{u) for varying shower age and particle energy is 
shown in fig. 5, and is found to be almost independent of these parameters. 
There is a slight broadening of the distribution at small shower ages, but 
this corresponds to the beginning of the shower when the number of charge 
particles is small. 




Fig. 5. Left: Variation of the distribution of logj^g ^ with shower age (for 1 TeV 
primaries). Right: Variation of the distribution of log^^Qii with particle energy. The 
distribution is sUghtly broader at small shower ages and at large energies, however 
the discrepancy is not very significant. Moreover, these regions of the parameters 
space correspond to a small number of charged particulars. 

The distribution of dlogiQ^u), assumed to be independent of shower age or 
primary particle energy, is shown in fig. 6 and can be modelled by the expres- 
sion 

^ (Po + Pi X logio U + P2X arctan(logio u - pa)) 



d loEi n u 



A X exp 



with: 



Po = l-55, pi=0.29, P2 = 2.5, = 0.73 for E< 300MeV 

Po = l-16, pi=0.50, P2 = l-8, P3 = 0.57 for E> 300MeV(12) 

The distribution dN/dlogiQ{u) for different primary energies and for charged 
particles energies E < 300 MeV (left) and E > 300 MeV (right) is shown in 
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fig. 6 with the analytical formula from eq. 11 superimposed (solid red line). A 
good agreement up to a few % is observed in the central part of the distribu- 
tion, encompassing the majority of particles in the shower. 




Fig. 6. Distribution of rescaled angle u for different primary energies ranging from 
10 GeV to 10 TeV, and for charged particle energies E < 300 MeV (left) and 
E > 300 MeV (right). 



1.3 Lateral distribution of particles in a shower 




Fig. 7. Definition of coordinate system attached to each electron/positron in the 
shower. 

The lateral distribution of the charged particles (electron or positron) in a 
shower are given in a system of coordinates attached to each charged particle 
in the shower as introduced in [10]. The X-direction is defined as the projection 
of the particle flight direction on the plane perpendicular to the shower axis, 
and the Y direction is perpendicular to it (and therefore perpendicular to 
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the particle velocity) (see fig. 7). The mean particle displacement (xg) with 
respect to the shower axis along the X axis is non-zero, whereas the average 
displacement along the Y axis (7/e) is zero for symmetry reasons. 

Assuming a factorisation of the expressions of mean displacement and spread 

in terms of energy, particle angle and shower age, the following expressions are 
found to properly describe the particle spread (expressed in units of g x cm~^): 

(xe) =24.3 X exp(0.63 x logw + 0.025 x log^ w) 

x(l + 0.2 X s^-^) X (1 - exp(-0.11 x (13) 
21 

X- X exp(0.47 X log^Mev - 0.023 x logE^.y) 

-C'MeV 

(7^^ =28.7 X exp(0.85 x logw) x (1 + 0.9 x w^-^) 

X (1 + 0.37 X s^-^) X (1 - exp(-0.03 x s°-^)) (14) 
21 

X- X exp(0.55 X log^Mev - 0.028 x log^J^ev) 

-C/MeV 

(?/e)=0 (15) 

cTj,, =2.65 X exp(0.03 x \ogw) x (1 + 0.51 x w°-^) 

X (1 + 0.2 X s*-^) X (1 - exp(-0.2 x s)) (16) 
21 

X- X exp(0.675 x log£;MeV - 0.035 x logEl^^^) 

-C/MeV 

The above expression is a good description of the mean lateral displacement 
and spread at the level of 10%. 

The reduced variables and Y^. are used for simplicity: 

Xr = Yr = ^ (17) 



were (xg), (t^^ and are obtained from eq. 13 to 16. The distributions of X^ 

and Yr- for 1 TeV gamma-ray showers are shown in fig. 8. The average values 
and RMS of X^ and Y^ are respectively and 1, as expected if equations 
13-16 are correct. As expected, the Y^ distribution is symmetric while the X^ 
distribution is not. 

The Yr distributions can be modelled by the following expression: 
dA^ / F2 \ 
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Xf Lateral Distribution 



Yr Laterai Distribution 




Fig. 8. Distributions of reduced lateral displacement Xr and Yj. (eq. 17) for 1 TeV 

showers, along X and Y axis. As expected, the distribution along Y axis is symmet- 
rical while the distribution along the X axis is not. 

A similar expression can be derived for Xr taking into account the non-centred 
position of the distribution maximum, and with different coefficients on both 
sides. For the sake of simphcity, a tabulated version of this distribution is used 
in the model production instead of a complicated analytical expression. 

The dependency of the reduced Xr and Yr lateral distributions with shower age 
and particle angle is shown in fig. 9. The Xr and Yr distributions remain stable 
enough to be considered as independent of these parameters in the shower 
modelling. The small variation seen in the distribution of Xr with particle 
angle affects mainly particles with a very large angle which are anyhow not 
likely to reach the telescopes. 

Figure 10 shows that the Xr and Yr distributions are mostly uncorrelated and 
can therefore be considered as independent. 



2 Model Generation 



The semi- analytical description of shower development described in the pre- 
vious section can be used to construct a shower model, i.e. a prediction of the 
light distribution on the ground for a given primary particle energy, direc- 
tion, impact parameter and development depth. This section describes how 
the shower image model is constructed from the distributions derived in the 
previous section. 
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Fig. 9. Evolution of the reduced lateral distributions {Xr- left, Y^rright) with shower 
age (top) and reduced particle angle (bottom). The strongest dependency is seen in 
the dependency of Xr on the reduced angle w and remains acceptable. 



2. 1 Principles 



The light density due to a shower in the camera can be calculated by an 
eight-fold integral: 

• integral over altitude z or depth t (longitudinal development of shower). 

• integral over energy of the electron/positron in shower. 

• integral over electron direction with respect to the telescope {u and 0). 

• integral over electron position with respect to its direction (X^ and Yr). 

• integral over Cherenkov photon wavelengths. 

• integral over Cherenkov photon azimuthal angle around the electron (the 
angle between the electron and the Cherenkov photon being fixed for a given 
electron energy). 
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Fig. 10. Bi-dimensional reduces lateral position distribution. There is a little corre- 
lation between the and the Y^. distribution. 

I{x,y) = JdzJdEx ^{t,E) x ^(y) 

J duxFME,s)) (19) 

J dXrJ dYrFxY{Xr,Yr,E,S,u) 
xCol{z, Xr, Yr, U, (j), (j)ph) 

Where: 

• dNe/dE{t, E) is the energy dependant longitudinal distribution of charged 
particles in the shower (from eq. 1 and 5) as function of atmospheric depth 
from first interaction 

• Fu{u{E,s)) is the normalised angular distribution of particles (from eq. 4 
and 11) 

• Exy{E, s, u) is the normalised lateral distribution of particles, (from eq. 13 
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to 16) 



• 1/A^ X (]?n^/ (cos dz dA) is the Cherenkov photon production rate (per unit 
of vertical track length and emitted photon wavelength) for an electron angle 
9 with respect to the shower axis 

• exp(— r(z, A)) is the atmospheric absorption 

• QeffW is the detector quantum efficiency (multiplied by mirror reflectivity 
and other wavelength-dependent transmission coefficients in the detector) 

• Col{z, X, Y, u, 0, (pph) is the average geometrical collection efficiency for pho- 
tons emitted by a electron at position (X, Y, z) with direction defined by 
(m, 0), and with an azimuthal photon angle (j)ph around the electron. 

In addition to the aforementioned ingredients, instrumental effects such as 
instrument point spread function and electronic response of the camera, in- 
cluding in particular trigger response and integration duration, have to be 
taken into account in the calculation. 

These effects, as well as the geometric light collection efficiency are obtained 
from a detailed simulation of the telescope response and parametrised in look- 
up tables. Atmospheric absorption of light, wavelength-dependent quantum 
efficiency and reflectivity used in the model generation are also implemented 
as look-up tables. 



Fig. 11. Definition of integration range for electron azimuthal angle (see text). 

In practice, for a given set of parameters, the relevant altitude range (in which 
emitted photons can reach the telescope) is flrst determined. Within this range, 
the shower development is divided into altitude slices (typically 20 slices). For 
each slice, the total number of charged particles in the shower is computed 
from eq. 1, and the particles are distributed in energy bands using eq. 5. For 
each energy band, the average energy, Cherenkov photon emission angle and 
emission yield (integrated over atmospheric transmission, quantum efficiency 
and dish reflectivity) is derived. The average shower is then further divided into 
flying particle angular bands and lateral displacement (eq. 7 to 18). For each 
considered position and direction, the geometric collection efficiency is then 



Telescope ring 
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taken into account to compute the contribution of this small shower subset. 
At each calculation step, a precise estimation of the required integration range 
is performed to avoid spending too much time in sampling parts of the shower 
whose Cherenkov hght do not reach the telescope mirror. 

An example for such a determination (using a geometrical construction) is 
shown in fig. 11: The projection plane (ground) is taken perpendicular to the 
shower axis, which is assumed to intersect the ground at point S. In order to 
avoid correlations with electron position and direction, the azimuthal range 
is determined in the shower direction frame (i.e. for a fixed shower and a 
telescope moving around the shower). The electron, at a fixed altitude z, with 
an energy E^, is assumed to be at a distance R with an angle 6 from the shower 
axis. Point E in the figure is the projection on the ground of the electron 
position. Its impact on the ground lies on the circle of radius Dg = zta.n9 
around the point E. Any Cherenkov photon emitted by this electron (with 
an angle 6c with respect to the electron) will fall within the ring defined by 
the radii Z}™" = ztan{9 — 9^) and D"'"-^ = 2;tan(^^ + 9^). Intersections of 
this ring with the possible telescope position (circle centred on S) given the 
allowed azimuthal range [0mm, 0moa;]- A similar approach is used to define the 
integration range for the other variables. 



2.2 Parameter space 

Models are generated for: 

• 40 different zenith angles 9z 

• 40 impact distances ranging from to (400 m)/ cos(^2) from telescope 

• 65 different energies from (50 GeV)/ cos(6'2) to (20 TeV)/ cos(^^) 

• 6 first interaction depths from to 5 Xq. 

The shower images are always generated on- axis (i.e. for a 7-ray at the centre 
of the camera). For a perfect telescope, a change of direction will result in a 
simple translation in the camera frame, that can be applied later, in the fit 
procedure. In a more realistic telescope, the broadening of the point spread 
function at large off-axis angles needs to be taken into account. This is cur- 
rently not needed for the H.E.S.S. telescopes, which use a Davis-Cotton optical 
design offering a degradation of the PSF from 0.25 mrad (RMS) at the centre 
to ~1 mrad at the edge of the field of view, however always smaller than the 
pixel sizc[12]. In total, 624000 templates are generated. The model generation 
procedure requires about 50 to 100 day-machine computing time on recent 
desktop computers and can be easily parallelled. The resulting shower images 
are stored in a ROOT binary filc[13] for later use. 
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2. 3 Example 



The output of the model generation procedure is a bank of two dimensional 
shower images in the frame of a perfect camera, with very small (0.01°) pixel 
size. For a given primary particle set of parameters, a predicted image is com- 
puted with an interpolation procedure on the generated images in a 4 dimen- 
sional space (energy, impact distance, primary interaction depth and zenith 
angle). Shower direction and azimuthal angle are then taken into account as a 
rotation and a translation in the camera frame to compute the final predicted 
image. Examples of such two dimensional shower images are shown in fig. 12 




Fig. 12. Model of a 1 TeV shower started at one radiation length and falling 20 m 
(top-left), 100 m (top-right) and 250 m (bottom- left) away from the telescope. X 
and Y axes are in units of degrees in the camera frame. Bottom-right: same as 
top-right but with a first interaction point located deeper in the atmosphere (at 3 
Xq). Note that the vertical scale (image amplitude) differs. 
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2-4 Comparison with simulation 



A comparison of the image shapes between simulation and model prediction 
is shown in fig. 13 for 1 TeV gamma-ray showers. The images were calculated 
for a H.E.S.S. camera, with pixels of 0.16° diameter. The image length and 
width were estimated using the standard Hillas parametrisation applied to the 
predicted images. In each plot, the average value of the simulation is drawn as a 
black histogram, with error bars indicating the shower-to-shower fiuctuations. 
The model prediction (image of average shower) is represented by a solid thick 
red line. 
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Fig. 13. Comparison between 1 TeV simulated shower images at zenith (black) and 
model prediction (red). Top left: Image amplitude as function of impact distance 
for primary interaction point of one radiation length. Top right: Image amplitude 
as function of primary interaction depth for impact distance p < 80 m (core of light 
pool). Bottom left: Image length (top) and width (bottom), in units of milliradians, 
estimated with the standard Hillas parametrisation technique, and as function of 
impact distance. Bottom right: Angular distance, in units of milliradians, between 
the image centre of gravity and primary direction as function of impact distance. In 
each plot, the simulation is drawn as a black histogram, with error bars indicating 
the shower-to-shower fluctuations. The model prediction is represented by a solid 
thick red line. 

The agreement between the model and the average values of the simulation is 
excellent up to core distance of about 300 m, where some trigger effects can 
explain the differences: at this distance, the total image amplitude does not 
exceed 100 photoelectrons, distributed over many pixels. Showers fluctuating 
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up to higher intensities are more hkely to trigger the telescope, thus resulting 
in a higher average image amplitude in the simulation compared to the model. 

The bars in the simulation histograms (in black) are due to shower-to-shower 
fluctuations, which are not taken into account in the model generation. At 
smafl core distance, and when taking into account the evolution with primary 
interaction depth, the shower intensity fluctuates by about 20% from shower 
to shower in addition to the fluctuation related to the depth of flrst interaction. 



2. 5 Conclusions 

The work presented in this section results in the generation of a shower image 
model, which is nothing more than an accurate prediction of the expected 
Cherenkov light distribution in the camera for a given set of primary particle 
parameters. The key ingredients in the model generation are the inclusion of 
depth of flrst interaction as shower parameter (as the main source of shower-to- 
shower fluctuations), and the precise description of longitudinal, lateral, and 
angular distributions of charged particles in the shower. The corresponding 
shower model is constructed once for all and can be applied to various zenith 
angles, telescope impact distance or off-axis angle. 

An alternate brute force approach would be the generation of the shower im- 
age model from massive simulations performed at every energy, zenith angle, 
impact distance, off-axis angle and primary interaction depth. The size of the 
parameter space and the need for a smooth model (in order to avoid local min- 
ima that would hamper the quality of the flt procedure) makes this approach 
less effective that the generation of an intrinsically smooth shower model. 



3 Fit procedure 

Once a shower model has been obtained, actual images on the camera can 
be compared to the ones predicted by the shower image model for a given 
set of primary parameters. A minimisation procedure is then used to obtain 
the most likely parameters of the incoming particle (energy, direction, impact, 
depth of the first interaction) under the hypothesis that the particle is a 7-ray. 
The minimisation procedure involves a precise comparison of the intensity in 
each pixel of the camera with the prediction from the model (interpolated 
between grid points to the actual set of parameters). In order to take into 
account the Poisson nature of the detected number of photons in each pixel, 
a log-likehhood approach has been chosen. 
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3.1 Pedestal and calibration 

In the absence of Cherenkov light, each pixel of the H.E.S.S. camera is illumi- 
nated by a significant amount of Night Sky Background (NSB) light, which 
largely dominates over the electronic noise and represents a single photoelec- 
tron rate ranging from 40 MHz off the galactic plane to ~ 300 MHz in the 
most luminous parts of the galactic plane. The pedestal is defined for each 
pixel as the charge distribution collected in the absence of Cherenkov signal. 
It is determined[14] for each pixel and for each observation run, using a clean- 
ing procedure that identifies the position of the shower image in the camera 
and rejects the corresponding pixels. The pedestal width is the combination 
of the electronic noise and the night sky background, the later being largely 
dominant. Due to varying atmospheric conditions, rotation of the sky on the 
camer and instrumental effects that depend on particular on temperature, 
the pedestal position and width can vary with time, with timescales of the 
order of a few minutes. The evolution of pedestal position and width with 
time for each pixel is recorded for every observation run and used during the 
reconstruction describe below. 

3.2 Pixel log-likelihood 

The shower image model provides a density of Cherenkov light in the focal 
plane. A convolution by the pixel size is done when loading the model to 
compute the expected signal fi in each pixel. For the sake of simplicity, the 
dependence of the PSF across the field of view is not taken into account in 
this calculation. The probability density (likelihood) to observe a signal of s 
(in units of photoelectron) in a pixel for an expectation value n is given by the 
convolution of the Poisson distribution of the photoelectron number n with 
the photomultiplier resolution. The latter is well represented by a Gaussian of 
width y^cTp + na^ where cXp is the width of the pedestal (charge distribution 
under pure noise, including night sky background) and the width of the 
single photoelectron peak (photomultiplier resolution): 

^ E ^ ^^^^ e.p [-^^^T^^j (21) 

In eq. 21, the actual pedestal width dp is different for each pixel, and depends 
in particular on the level of night sky background (NSB) seen by the pixel. 
The photoelectron resolution a^, also specific to each pixel, is measured using 
calibration runs where the camera is illuminated with a low intensity flashing 

^ The H.E.S.S. telescope use a Alt-Az mount 
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LED. The use of measured values of cxp and o"^ is an important aspect that will 
explain the stability of the Model Analysis for varying level of NSB (Section 
4.12). 




Fig. 14. Probability to observe a signal of s (in units of photoelectron) in a pixel 
for an expected average value n, as function of s and fi, and for a pedestal width 
of 0.3 p.e. and a single photoelectron resolution of 40%, from equation 21. For low 
predicted or detected intensity, the simple Gaussian assumption (as used in a 
test) is not valid anymore. 

The shape of function 21 is shown in fig. 14. In order to obtain a variable that 
behaves asymptotically like a x^? we define the pixel log-likelihood 

In L = -2 X In P{s\fi,ap,a^) (22) 
3.3 Pixel log-likelihood expectation values 

In order to compare the observed signal s with a model prediction /i, the ex- 
pectation value of the log-likelihood under different realisations of the same 
shower image {i.e. being fixed and s fluctuating due to Poisson noise, elec- 
tronic and NSB fluctuations) needs to be computed. If the observed signal is 
only due to noise (night sky background), the probability density functions 
simplifies to a Gaussian of width ap. In this case, the average value of the pixel 
log-likelihood reads 

(In L) 1^=0 = ds \n L(s|/i = 0, ap) x P(s|/i = 0, ap) 
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= l + ln(27r)+ln(7j 



(23) 



In a similar manner, one obtains 

[in^ l)=(i + ln(27r) + In % 2 ^ (In L) = 2 (24) 



At high fi, the Poisson distribution can be replaced by a Gaussian of width 
y/]!, and the probability density function can be well approximated by the 
convolution of two Gaussian: 



P{s\ijL :$> 0, (7p, (T-y) pa . = exp 



/27r (a2 + /^(l + a2)) V 2 (^ci^ + ;,(1 + 

In this limit, the log-hkelihood behaves hke a x^: 

(In L) 1^ = 1 + ln(27r) + In (aj + //(I + a^)) , a^ln L) = 2 (25) 

This expression reduces to eq. 23 when /i — 0. Figure 15 shows the comparison 
of this analytical expression with a numerical integration of the average log- 
likelihood. The analytical expression is valid for fi = and, as expected, is 
also asymptotically valid for 3> 1. In the transition regime, the analytical 
expression slightly overestimates the likelihood value. 

In order to properly calibrate the average log-likelihood, the discrepancy be- 
tween the analytical expression and the numerical integration is computed for 
every expected amplitude /j, — and pedestal width cr^ (using a Monte Garlo 
simulation) and stored in a look-up table. This look-up table will be used at 
the end of the fit procedure (see below) to produce calibrated discriminating 
variables. 

3.4 Telescope log-likelihood 



Pixels belonging to a camera are assumed to be independent. We define the 
telescope log-likelihood as the sum over all pixels of the pixel log-likelihood: 

lnLtei= ln^i= H -2 X In P{si\i^i,ap,a^). (26) 

pixel i pixel i 
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Fig. 15. Average value and sigma of pixel log-likelihood as function of expected 
amplitude fi (see text) for different pedestal width (up = 0.5, l,and 2 p.e.). The 
solid line are the results of a numerical integration, whereas the dashed line are 
the analytical expression 25. Inset: difference between analytical expression and 
numerical integration. 



3.5 Reconstruction - Fit algorithm 



The model photon reconstruction relies on the pixel-per-pixel comparison of 
the actual shower images with the ones that are predicted for a given set of 
parameters. A minimisation procedure is used to find the best parameters 
(direction, impact, depth of the first interaction and energy). In contrast to 
Hillas parameters based reconstruction techniques, the raw image amplitudes 
are directly used, without any image cleaning procedure. All pixels are used 
in the fit, not only those close to the actual image. 

A Levenberg-Marquardt fit algorithm [15,16] is used to minimise the telescope 
log-likelihood (eq. 26). This algorithm is very efficient in the case the first and 
second derivative of the minimised function (log-likelihood or x^) can be ex- 
pressed analytically and depend mostly on the first derivative of the model (the 
second derivative being negligible). This is in general valid when the minimised 
function is a quadratic form [x^ or similar) and when the model function ex- 
hibits a smooth behaviour when varying the model parameters. The algorithm 
assumes a quadratic form and uses the inverse matrix of second derivatives 
to estimate the position of the minimum. In order to improve convergence 
stability, a combination of pure steepest descent and parabolic assumption is 
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used with weights that vary during the minimisation, depending on whether 
the quadratic assumption is vahd or not. The Levenberg-Marquardt fit algo- 
rithm always converges - sometimes to a local minimum - and provides reliable 
estimates of the model parameters uncertainties, taken from the correlation 
matrix. 

An important issue in the fit algorithm is the starting point. The standard 
Hillas reconstruction technique[17] with different sets of image cleaning pa- 
rameters is used to derive a handful of possible estimates. The estimate that 
provides the best initial log-likelihood is chosen as starting point of the min- 
imisation. 

In our case, about 40 iterations arc needed for the algorithm to converge, 
whereas simpler algorithms such as Minuit[18] needed in average about ten 
times more. On recent desktop computers, the reconstruction of a single event 
takes about 50 to 100 ms. 

The actual telescope refiectivity, as measured from ring-shape images of muons 
passing through the telescope[14], is used to scale the model prediction to 
the observation conditions, on a run-by-run basis. Since the optical efficiency 
is taken into account directly at the reconstruction level, no further energy 
correction is needed (which is completely different to the way the optical 
efficiency is handled in Hillas parameters based analyses). 

The output of the minimisation procedure are: 

• Best guess of the 6 shower parameters: direction (2 parameters), impact (2 
parameters), depth of the first interaction and energy 

• Correlation matrix and therefore uncertainty on these best fit parameters 

• Final log-likehhood value 



3. 6 Goodness of fit 

Since Atmospheric Cherenkov Telescopes are background dominated systems, 
the performance of any analysis is mainly driven by its ability to discriminate 
between gamma-ray induced showers and the much more numerous hadronic 
background. In the Model Analysis, a goodness- of- fit approach is chosen to 
compare the model prediction and the actual shower images, in order to check 
the compatibility of the recorded events with a pure 7-ray hypothesis. The 
goodness- of- fit is defined as a normalised sum over all pixels of the difference 
between the actual pixel log-likelihood and its expectation value, properly 
normalised: 
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[in L(si\ Hi) - {In L) 



pixel i 



(27) 



V2 X NdF 



where NdF is the number of degrees of freedom (number of pixels - 6). The 
goodness of fit behaves asymptotically like a and provides therefore a mea- 
sure of the fit quality. This will be used later on the hadronic discrimination. 
By construction, if the pixels likelihood behave like independent random vari- 
ables, G is expected to behave like a normal variable: 



3. 7 Goodness of fit calibration 

The goodness of fit average value rehes on the assumption of a Gaussian 
pedestal (of width Up). In the absence of night sky background, the pedestal 
reduces to the electronic noise and the assumption is always valid[14]. The 
H.E.S.S. camera use two different gains: a high gain, with dynamical range 
from to 200 photoelectrons with single photoelectron resolution and an elec- 
tronic noise representing about 0.2 p.e., and a low gain, with dynamical range 
up to 2000 photoelectrons but worse resolution. In the low gain, where the 
electronic noise represents about 2.5 p.e., the pedestal is also Gaussian to a 
good approximation. 

In the high gain channel, figure 16, reproduced from [14], shows that the Gaus- 
sian assumption is also almost valid for high night sky background (above 
~ 150 MHz), but not in the intermediate regime where the pedestal shape ex- 
hibits two peaks, similar to a single photoelectron spectrum, and produced by 
fractions of single photoelectron pulses falling by chance within the acquisition 
window. 

Simulation of 7-ray spectra were conducted with different levels of night sky 
background (from 10 to 800 MHz). Figure 17 shows the average goodness of fit, 
computed on a single telescope basis (to keep a constant number of pixels), 
and as function of the night sky background level for the two gains of the 
H.E.S.S. camera. 

The goodness of fit calculated with low gain channels on each pixel appears 
to be rather stable with night sky background, thus confirming the validity 

of the Gaussian assumption for this channel. In contrast, the goodness of fit 
calculated with high gain channels exhibit a strong deviation at low night sky 
background level, where the pedestal is not Gaussian anymore. Simulations 
were used to derive a correction table as function of night sky background 



(G) = and a^G) = 1 



(28) 
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Fig. 16. Evolution of pedestal shape with increasing night sky background level 
(solid line) and Gaussian fits to the distributions (dashed lines), using only the high 
gain channel. From [14]. 
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Fig. 17. Average value of goodness-of-fit per telescope as function of night sky 
background level, and as obtained from simulations, before (left) and after (right) 
goodness calibration. 



level in each pixel independently. This table is not used during the reconstruc- 
tion, and is only applied to the final goodness of fit estimator. The resulting 
goodness-of-fit is shown in fig. 17, before correction (left) and after correction 
(right). After calibration, the goodness-of-fit is stable for a night sky back- 
ground variation from 10 to 200 MHz, which corresponds the bulk of the 
H.E.S.S. observations. 
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3.8 ShowerGoodness and BackgroundGoodness 

The discrimination between 7-ray induced showers and hadron induced ones 
can use several distinct features (fig. 18): 

• Hadron induced showers are more irregular, and contain several electromag- 
netic sub-showers initiated in particular by disintegration of neutral pions. 
As a consequence, the image in a Cherenkov camera often exhibits several 
clusters separated apart (fig. 18, left) 

• In addition, the hadronic component of the shower itself emits a low in- 
tensity Cherenkov light spread over a large fraction of the camera. This 
emission, denoted as Hadronic rain, is in general eliminated by cleaning 
procedures used in standard reconstruction techniques, (fig. 18, middle). 

• Finally, the charged nucleus entering the atmosphere can emit a Cherenkov 
light before interacting. This emission is produced very high in the atmo- 
sphere and therefore it generates, in the camera, a faint spot close to the 
shower direction. 




Fig. 18. Examples of image topologies that can discriminate between a 7-ray induced 
showers and an hadronic one. Left: image subdivided in several clusters correspond- 
ing to electromagnetic sub-showers. Middle: hadronic rain emitted by the hadronic 
component of the shower. Right: Cherenkov emission from the primary particle (nu- 
cleus) entering the atmosphere, before the actual shower development. The red cross 
denotes the centre of the camera. 

In order to fully exploit the differences between the 7 and hadron induced 
showers, individual pixels contributing to the goodness-of-fit (eq. 27) are clas- 
sified into two different groups at the end of the fit: 

• Pixels belonging to the shower core, defined as pixels whose predicted ampli- 
tude is above 0.01 p.e., are grouped together with three rows of neighbours 
around them to construct a variable named ShowerGoodness (SG) in a sim- 
ilar way to eq. 27. These pixels are selected at the end of the fit procedure 
to avoid changes of the number of degrees of freedom during the fit. Due to 
the large reduction of the number of degrees freedom, this variable is more 
sensitive than the Goodness to discrepancies between the model prediction 
and the actual shower images. 
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• Remaining pixels, denoted as background pixels, arc grouped together to 
construct a variable named BackgroundGoodness (BG), which is sensitive 
to hadronic clusters outside the main image, hadronic rains and other ir- 
regularities. 



4 Results 

4-1 Data sets 

For a detailed comparison of the Model Analysis with previously existing 
reconstruction techniques, two data sets are used: 

• A data set {PKS Flare) of 7.7 hours (live time) of data taken on the blazar 
PKS 2155-304 during flaring period in July 2006, with 4 telescopes, yield- 
ing a test sample of more than 10 000 7-rays with very small background 
contamination, at relatively low zenith angles (typically 20°). The average 
night sky background for this data set corresponds to a single pixel trigger- 
ing rate of 40 MHz, representative for extragalactic space. Since the data 
has been taken within 3 days, the system state can be considered as stable. 

• A large set ( Crab Full) of ~ 34 hours (Uve time) of data taken on the Crab 
Nebula, from 2004 to 2008, with 4 telescopes and larger zenith angles (from 
40° to 60°). This data set yields a sample of ~ 10 000 excess events. The 
average NSB for this data set, representative for galactic plane, is more than 
twice higher (~ 100 MHz) than for extragalactic space. 

The model analysis will be compared to the standard Hillas parameters based 
reconstruction in use in the H.E.S.S. collaboration[17]. Two sets of cuts will 
be defined: Hillas 60 will denote a minimum image amplitude of 60 photo- 
electrons (p.e.), same as for the Model Analysis, and Hillas 200 will use a 
larger minimum image amplitude of 200 p.e.. When needed, Hillas Std and 
Hillas Hard will denote configurations used in the Crab publication [17] with 
respectively 80 and 200 p.e. minimal image amplitude. 

4-2 Shape cuts 

In order to have well reconstructed showers, the following shape cuts are used 
throughout the current section (unless specified differently): 

• A minimum image amplitude of 60 photoelectrons per telescope. 

• A maximal nominal distance (distance of the shower image centre to the 
centre of the camera) of 2° for a camera radius of 2.5°. The nominal distance 
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cut removes truncated images, close to the camera edge, that often lead to 
misreconstruction of shower direction. 
• At least images from two different telescopes passing the previous cuts to 
ensure stereoscopic reconstruction. 

These cuts will be used both for the Model Analysis and for the Hillas param- 
eters based reconstruction techniques. The later needs an additional cleaning 
procedure. We use a two-level filter cleaning[17] with pixel thresholds of re- 
spectively 5 and 10 p.e.. The image amplitude is computed after cleaning. The 
Model Analysis doesn't require any cleaning procedure, so the image ampli- 
tude for the model photon reconstruction is slightly larger than for the Hillas 
analysis (as pixels belonging to the shower image tail are taken into account 
for the model reconstruction and not for the Hillas reconstruction). 

4-3 Hadron discrimination 




Shower Goodness Background Elfioiency 



Fig. 19. Left: Distribution of Shower Goodness for real data taken on the blazar 
PKS 2155-304 (blue point for excess events, grey triangles for background events), 
compared with a simulation (red histogram) with a similar night sky background 
levek Right: 7-ray efficiency as function of background rejection for a selection 
based on Shower Goodness, compared with selection achieved in a standard Hillas 
parameters reconstruction at the same minimum image amphtude (60 p.e.) and 
with the Mean Scaled Width and Mean Scaled Length variables. 

The distribution of ShowerGoodness (SG) in the shape cuts of section 4.2 is 
shown in fig. 19 for the data taken on the blazar PKS 2155-304, and for a 
simulation with 40 MHz of night sky background noise, corresponding to the 
average night sky luminosity around PKS 2155-304. The Monte Carlo dis- 
tribution (red line) and the distribution for the PKS 2155-304 excess events 
(blue filled circles) are in overall good agreement, and are very different from 
the distribution obtained for background events (grey triangles), thus con- 
firming the discrimination capabilities of the ShowerGoodness variable. The 
small shift between the two distribution is responsible for a systematic error 
in acceptance determination of 4% for a cut at SG < 0.6. 
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Fig. 19, right, shows the 7-ray efficiency as function of background rejection 
for a selection based on ShowerGoodness onl}{f], after shape cuts, compared 
to what is achieved in the standard Hillas parameters based reconstruction 
(at the same minimal image amplitude, and using the Mean Scaled Width 
as varying selection variable). A large squared angular distance cut of 6'^ < 
0.04 deg^ was used to ensure that all events are accepted (extended source 
assumption) and to decorrelate angular resolution performance(sec. 4.8) from 
hadronic discrimination. The model reconstruction provides, for a given 7-ray 
efficiency, a much better background rejection than Hillas parameters based 
reconstruction. Moreover, as will be shown in the following sections, the Model 
Analysis provides additional discriminating variables that improve further the 
sensitivity. 

In the standard configuration of the Model Analysis, a cut ShowerGoodness 
SG < 0.6 will be used as main discriminating parameter. This cut retains 70% 
of 7-rays and rejects more than 95% of background events, yielding a quality 
factor 

Q= ^4 (29) 

V ^hadrons 



4.4 Primary interaction reconstruction 



The depth of the ffist interaction being a parameter of the model, is also a 
direct product of the reconstruction procedure (instead of being calculated, 
as in the Hillas parameters based reconstruction, from the shower maximum). 
Fig. 20 shows the ability of the Model Analysis to reconstruct the depth of the 
first interaction with almost no bias and a resolution of 0.7 Xq (slightly worse 
at large zenith angles). This is better than Hillas parameters based analyses, 
which obtain resolution not better than 1 Xq in the best cases. 

A comparison of the actual depth of the ffist interaction obtained with real 
data and with simulation is shown in fig. 21, in shape cuts and SG < 0.6, 
under a point-like source assumption (squared angular distance cut of 9'^ < 
0.01 deg^). The agreement between data and simulations is excellent. The 
PKS 2155-304 data are compatible with a 0.7 Xq resolution, and the Crab data 
with a slightly worse (0.9 Xq) resolution due to a larger average zenith angle. 
Moreover, since the distribution for OFF data is significantly different, the 
depth of the ffist interaction can be used to improve the 7-hadron separation. 
The irregularities in the OFF distributions are artifacts of the reconstruction 
procedure: models are only available for depths of 0,1,2,3,4 and 5 Xq and 



^ several other variables, described later in this paper, can provide additional re- 
jection 
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Fig. 20. Reconstruction of the depth of the first interaction for a 7-ray spectrum of 
differential index 2.0 at zenith. Left: Reconstructed depth of the first interaction (in 
radiation length) versus true one. Right: profile showing the linear relation between 
true and reconstructed depth of the first interaction. 




Fig. 21. Distribution of depth of the first interaction for real 7-ray (from observations 
on PKS 2155-304 for the left figure, and on the Crab nebula on the right) compared 
to Monte-Carlo simulations, respectively at zenith with differential index —3.4 and 
at a zenith angle of 46° with differential index of —2.6. The agreement between the 
simulation and the actual data is very good. 

are linearly interpolated between these values (and extrapolated above 5 Xq). 
Poorly constrained showers tend to accumulate at the grid model points and 
avoid interpolated values. 



In order to improve 7-background separation, a cut in reconstructed primary 
interaction depth —IXq < tQ < AXq is used. This cut retains 90% of 7-rays 
(resp.91%) and rejects 30% of hadrons (resp. 25%) , yielding a quality factor 
Q = 1.07 (resp. 1.06) respectively for the Crab and PKS 2155-304 data sets, 
leading to an improvement of the signal over background ratio by about 20%. 

At low energies, the electron background becomes dominant for atmospheric 
Cherenkov telescopes. It is usually considered as irreducible, since electron 
induced showers are almost identical to 7-ray induced one. However, electron 
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Fig. 22. Distribution of the depth of the first interaction for a 7-ray spectrum of 
differential index 2.6 at zenith (in red) compared to electrons with same spectrum 
(in blue): electron induced showers start to emit Cherenkov light a little bit earlier 
in the atmosphere. 

induced showers start to produce Cherenkov light a little bit earlier in the 
atmosphere (about a radiation length before 7). This is confirmed by simu- 
lation shown in fig. 22: the distribution of reconstructed primary depth for a 
simulated electrons spectrum shows a peaks centred around Xq, earlier than 
for 7's. The difference is not big enough to distinguish between electrons and 
gammas on an event by event basis, but a statistical discrimination seems to 
be possible. 

4-5 Analysis Configurations 

The 7-hadron separation is essentially based on the Shower Goodness and pri- 
mary interaction depth variables. For completeness, two additional cuts are 
used: a cut on BackgroundGoodness {BG < 2), with 7-ray efficiency close 
to 100%, is designed to remove a small number of shower misreconstructed 
at very large distance from the camera centre (4 degrees or more). A second 
cut on the Goodness (G) variable (eq. 27), in which the model is replaced 
by a null (yU = 0) assumption, removes images that are consistent with pure 
night sky background noise, and reduces the effect of night sky background 
inhomogeneities across the field of view. 

The standard cuts for the Model Analysis are defined as: 

• A minimum image amplitude of 60 photoelectrons per telescope 

• A maximal nominal distance (distance of the shower image centre to the 
centre of the camera) not more than 2° 

• At least two telescopes passing the previous shape cuts 

• A maximum ShowerGoodness (SG) of 0.6 

• A reconstructed primary interaction depth to between -1 and 4 Xq. 
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• For a point-like source, a squared angular distance cut of 9'^ < 0.01 deg^, 
independent of zenith anglq£j 

Two additional sets of cuts are defined as reference point for future H.E.S.S. 
publications: the Faint source configuration is optimised for source fainter than 
a few percent of the Crab flux. The Loose cuts configuration is designed to 
maximise the 7-ray efficiency for strong sources at the expense of a poorer 
background rejection. For analysis of extended source, the 9^ cut is usually 
replaced by a selection that encompasses the whole source. Cuts for the three 
aforementioned configurations are summarised in Tab. 1. of cuts 



Name 


Min. 


Max. Nom. 


#Tels 


SGmSuX 


to 


n2 

'^max 




Charge 


Distance 












(p.e.) 


(deg.) 






(^0) 


(deg2) 


Standard 


60 


2 


2 


0.6 


[-1,4] 


0.01 


Faint Source 


120 


2 


2 


0.4 


[-1,4] 


0.005 


Loose Cuts 


40 


2 


2 


0.9 


N/A 


0.0125 



Table 1 

Reconstruction configurations: Minimum charge, maximum nominal distance, max- 
imum ShowerGoodness, primary interaction depth range and maximum squared 
angular distance from the reconstructed shower position to the source. A minimum 
of two telescopes passing the per-telescope cuts, on image amplitude and distance 
from the centre of the field of view, are also required. 

4-6 Effective Areas 

The effective function of energy for the Model Analysis at zenith is 

shown in fig. 23, compared to areas obtained with Hillas parameters based 
reconstructions. The reconstruction efficiency for the standard configuration 
is similar to the values obtained with Hillas reconstruction with a minimum 
image amplitude of 60 p.e.. As expected, the loose configuration has a larger 
effective area and a lower threshold. In most of the H.E.S.S. publications so 
far, a minimum image amplitude of 200 p.e. was used, yielding a much larger 
threshold of ~ 300 GeV at zenith. Models are currently generated only up 
to 20 TeV, thus leading to a loss of acceptance at high energy. In addition, 
showers above 10 TeV can reach the ground, resulting into large fiuctuations 
that are not fully reproduced by the model. 



Since the angular resolution degrades at large zenith angle, a varying cut could 
perform better. This is not specific to the Model Analysis and applies to other 
reconstruction techniques as well. 
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Fig. 23. Effective area as function of energy, at zenith, compared witli the values 
obtained for the standard and hard cuts Hihas parameters based analyses. 

Performances obtained on the Crab Nebula and on the Blazar PKS 2155-304 
are shown in tab. 2. The Model Analysis yields a 7-ray efficiency similar to 
the 60 p.e. Hillas reconstruction, but with a background rejection improved 
by a factor of 6.2 (PKS 2155-304) to 6.8 (Crab Nebular), yielding a signal 
over background ratio similar or better to the one obtained with hard cuts 
Hillas reconstruction. As a result, the significance obtained with the Model 
reconstruction is larger, and the sensitivity improved by a factor of more than 
2. The faint source configuration provides an additional factor of ~ 2 in S/B 
compared to standard configuration. 



Ji-.l Energy resolution 



The energy resolution (in cuts) as function of energy is shown in fig. 24 for 
zenith, where the energy resolution is defined as the RMS of the AE /E distri- 
bution. The energy resolution is better than 15% for the whole energy range 
(from 80 GeV up to more than 20 TeV), with biases not exceeding 5% in the 
central range. The very central energy range (500 GeV to more than 10 TeV), 
the energy resolution is better than 10% and reaches values as low as 7 to 8%. 
Larger energy biases appear at very low energy (up to 20% at 80 GeV), due to 
trigger selection effects. These biases are however smaller than those obtained 
with an Hillas parameters based reconstruction. In a similar manner, negative 
energy biases appear at very high energies, due to very distant high energy 
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Table 2 

Number of excess events and significances obtained with the Model Analysis on 
some H.E.S.S. sources compared to standard reconstruction techniques. The pub- 
lished results on the Crab Nebula[17] where obtained with a minimum image am- 
plitude of 80 photoelectrons, yielding a better S/B ratio than results obtained with 
a 60 pe cut, although with similar significance. 

showers reconstructed closer to the telescope£l] and with an underestimated 
energy. In the medium energy range (500 GeV to a few 10 TeV), the Model 
Analysis largely outperforms standard reconstruction techniques. 



4-8 Angular resolution 



The angular resolution, defined as the 68% containment radius, is shown in 
fig. 25 for showers at zenith (left), and is of the order of 0.06° in most of the 
energy range, which is much better than values obtained with Hillas param- 
eters based analyses (using algorithm 1 from [19]) for similar minimal image 
amplitude. An immediate consequence is an improved sensitivity for point like 
sources and improved morphological studies of extended source. The angular 
resolution is stable for zenith angles up to 50° (fig. 25, right), ans rises slowly 
up to very large zenith angles. The degradation observed for very large zenith 
angles is much larger for Hillas parameters based reconstruction techniques. 

The superior angular resolution of the Model Analysis is demonstrated in 
fig. 26 on data taken on PKS 2155-304. The 6"^ distribution is twice as more 

Very distant shower produce almost parallel images in the camera, which intro- 
duces a degeneracy in the reconstruction. 
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Fig. 24. Energy resolution (main plot) and bias (inset) as function of energy, at 
zenith, compared with the values obtained for the standard and hard cuts Hillas 
parameters based analyses (resp. blue circles and triangles). 
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Fig. 25. Angular resolution (defined as the 68% containment radius) as function of 
energy, at zenith, compared with the values obtained for the standard and hard cuts 
Hillas parameters based reconstruction techniques. Right: Average angular resolu- 
tion for a spectrum, as function of zenith angle. 

peaked for the Model Analysis (right), compared to the Hillas reconstruction 
(left). The same behaviour is observed at larger zenith angles for the Crab 
Nebula (fig. 27). 

Detailed numerical comparison shown in tab. 3 demonstrates that, on average, 
the Model Analysis angular resolution performances with a minimal image 
amplitude of 60 p.e. are equally good as the Hillas Analysis with a much 
higher minimal image amplitude of 200 p.e.. 
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Fig. 26. Squared angular distribution [6'^) obtained on the blazar PKS 2155-304 with 
the standard Hillas parameters based analysis with a minimal image amplitude of 
60 p.e. (left) and with the Model Analysis (right). 




Fig. 27. Squared angular distribution (0^) obtained on the Crab Nebula with the 
standard Hillas parameters based analysis with a minimal image amplitude of 60 
p.e. (left) and with the Model Analysis (right). 
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Table 3 

Angular resolutions obtained on real data (Crab Nebula and PKS 2155-304 for the 
Model Analysis and Hillas parameters based analyses. 
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4-9 Uncertainty on parameters 



The Model Analysis provides uncertainty estimation for each reconstructed 
parameter, as byproducts of the correlation matrix. This can be used to select 
in a simple and natural way a sub-sample of the events with improved angular 
or energy resolution (for more precise morphological or spectral analysis). 

An example of such a subset selection is shown in fig. 28, where an additional 
selection on the energy uncertainty dlnE < 0.04 is applied to the data. The 
energy resolution becomes better than 8% from 100 GeV up to 10 TeV with 
values as low as ~ 6 % at a few TeV. The prices to pay for such a golden 
events selection is an increase of energy threshold and a reduction of statistics 
(fig. 28, right) by a factor of ~ 2 above 1 TeV. 




Fig. 28. Effect of an additional selection on the energy uncertainty dlnE < 0.04 on 
the energy resolution (left) and effective area (right). 

In a similar manner, fig. 29 shows the effect of an additional selection on the 
direction uncertainty dDir < 0.03° on the angular resolution. It is possible to 
achieve a resolution as good as 0.03° in the TeV energy range at the expense 
of a factor less than 2 in statistics. This can greatly improve the ability to pin- 
point the source position in challenging regions such as the Galactic Centre. 
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Fig. 29. Effect of an additional selection on the direction uncertainty dDir < 0.03° 
on the angular resolution (left) and effective area (right). 
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This is demonstrated in fig. 30, where the same selection dDir < 0.03° is 
apphed to data from PKS 2155-304. The net result of the selection is an 
unprecedented resolution ergs = 0.055 and ergs = 0.1 with number of excess 
event being still around 10000 7. The background is also much more suppressed 
than the signal, resulting into a. S/B ratio of 454 compared to 236 before 
selection. The direction uncertainty could therefore serve as an additional 
discriminating parameter, but its effect would not be independent of energy. 
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Fig. 30. Effect of an additional selection on the direction uncertainty dDir < 0.03° 
on the squared angular resolution of events from the blazar PKS 2155-304. 



4-10 Combination of reconstruction methods 



In the early stages of the Model Analysis, it was realised that the Good- 
ness variable was almost uncorrelated with the Mean Scaled Width and Mean 
Scaled Length variables obtained with standard reconstruction techniques [20, 21] 
As a consequence, the combination of several reconstruction techniques was 
improving rejection capabilities and thus final sensitivity. 




Mean Scaled Width Mean Scaled Length 



Fig. 31. Distribution of Mean Scaled Width (left) and Mean Scaled Length (right) 
on real data taken on PKS 2155-304, for excess events (red), background events 
before selection on ShowerGoodness (blue) and after (black). All histograms are 
normalised to the same number of events. 

Fig 31 shows the behaviour of the Mean Scaled Width and Mean Scaled Length 
variable (from Hillas reconstruction) when imposing the ShowerGoodness cut. 
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The background data histograms (in blue) are much more extended than the 
excess events {gamma-rajs) histograms (red), thus confirming the rejection 
capabihties of the Mean Scaled Width and Mean Scaled Length variables. 
However, after imposing a selection on the ShowerGoodness {SG < 0.6), most 
of the discrimination capability vanishes (black histograms) : the Mean Scaled 
Width still provides only a very small additional discrimination, whereas the 
Mean Scaled Length does not anymore. The games of combining several anal- 
yses appears therefore not as promising as it was, when using a older version 
of the Model Analyis with less performance. 



4.11 Sensitivity to stars and broken pixels 
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Significance Distribution - Hillas 60 
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Significance Distribution - Model 
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Fig. 32. Significance Maps (left) and distributions (rigfit) obtained respectively with 
the Hillas Analysis (top) and Model Analysis (bottom), using Template Model back- 
ground subtraction. The Hillas Analysis is sensitive to a bright star ((^-Tauri, visual 
magnitude 2.97) in the field of view that causes a deficit in hadron acceptance. 

During H.E.S.S. observations, camera pixels illuminated by a bright star are 
turned off to avoid damage to the photocathode. In addition, some (in gen- 
eral less than 5%) pixels are removed from the analysis due to instrumental 
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malfunction (non working high voltage, noisy pixel, badly behaving analogue 
memory, ...). Missing pixels affect more Hillas parameters based analyses (as 
they produce truncated images) than the Model Analysis, which just ignores 
the missing pixels. 

The effect of bright stars is demonstrated in fig. 32, where the large Crab 
Nebula data set (32 live hours) was used to produce a significance map using 
the template model[22]. A bright star, ^-Tauri of visual magnitude 2.97, is 
present in the field of view. On average, about three to four pixels surrounding 
the star are turned off during observations. An artificial excess of events, up to 
the level of 8a, is visible in the Hillas Analysis significance map. This excess 
is due to a deficit in background events which translates into an artificial 
excess of events. The Model Analysis significance map shows no deviation at 
the same position, demonstrating that it is less sensitive to missing pixels. 
The artifact in the Hillas Analysis significance map weakens at larger minimal 
image amplitude (200 p.e.) as the shower images in the camera become larger 
and are therefore less affected by a few missing pixels. 



4-12 Sensitivity to varying night sky background 




Fig. 33. Left: Evolution of Shower Goodness distribution for simulated 7-rays with 
varying night sky background level. Right: Evolution of the mean Shower Goodness 
value with increasing night sky background level. 



The evolution of Shower Goodness distribution for simulated 7-rays (with 
photon index 2.0 at zenith) is shown in fig. 33. No strong evolution is seen 
up to 300 MHz. The bulk of the HESS observations correspond to an average 
NBS level of 100 MHz, varying between 40 MHz to 300 MHz in the most 
luminous parts of the Galactic Plane. 
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Conclusion 

We have presented a sophisticated 7-ray reconstruction technique for Atmo- 
spheric Cherenkov Telescopes, based on an accurate pixel per pixel comparison 
of observed intensity with a pre-calculated model. This significant improve- 
ment over earlier works leads to an improved angular and energy resolution, 
and background rejection, resulting in an increase of the sensitivity of the 
H.E.S.S. telescope system by a factor close to 2. The model analysis is less 
sensitive than reconstruction techniques based on Hillas parameters to in- 
strumental or environmental effects (missing pixels, night sky background 
variation across the field of view, . . . ), thus allowing an optimal use of the 
full dynamical range of the instruments. This novel reconstruction technique 
should also benefit to upcoming projects such as CTA and AGIS. 
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